/*
Copyright (c) [2019年5月1日] [吴超]
[MBT_Studio] is licensed under Mulan PSL v2.
You can use this software according to the terms and conditions of the Mulan PSL v2.
You may obtain a copy of Mulan PSL v2 at:
		 http://license.coscl.org.cn/MulanPSL2
THIS SOFTWARE IS PROVIDED ON AN "AS IS" BASIS, WITHOUT WARRANTIES OF ANY KIND,
EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO NON-INFRINGEMENT,
MERCHANTABILITY OR FIT FOR A PARTICULAR PURPOSE.
See the Mulan PSL v2 for more details.
*/
#pragma once

#include <向量.h>
#include <几何数据.h>
#include <矩阵.h>

#ifdef DEF_CUDA
	#include <SIMD/GPU_向量计算.cuh>
#else
	#include <SIMD/向量运算.h>
	#include <SIMD/CPU_向量计算.h>
#endif




#include <A_CUDA数学Core/aabbtree.h>
#include <A_CUDA数学Core/voxelize.h>
#include <A_CUDA数学Core/sdf.h>
#include <B_几何图形/convex.h>
//#include <B_几何图形/三角形.h>







typedef struct {
	vec3 pointA;
	vec3 pointB;
	vec3 pointC;
}S_三角形;





CUDA_CALLABLE inline float SampleSDF(const float* sdf, int dim, int x, int y, int z) {
	assert(x < dim && x >= 0);
	assert(y < dim && y >= 0);
	assert(z < dim && z >= 0);

	return sdf[z * dim * dim + y * dim + x];
}


CUDA_CALLABLE inline float f_SampleSDF(const float* sdf, const uvec3& dim, int32 x, int32 y, int32 z) {
	assert(uint32(x) < dim.x && x >= 0);
	assert(uint32(y) < dim.y && y >= 0);
	assert(uint32(z) < dim.z && z >= 0);

	return sdf[z * dim.x * dim.y + y * dim.x + x];
}







CUDA_CALLABLE Inline vec3 vec_Reflect(const vec3 vecI, const vec3 vecN) {
	//vecR = vecI - 2(vecI*vecN)vecN;
	//v' = v - 2 * (v*N) *N
	return vecI - (vecN * (vec_dot(vecI, vecN) * 2.0f));
}








template<typename T>
std::vector<T> f_Range(int64 s, int64 e) {
	int64 size = std::abs(e - s);
	std::vector<T> range(size);
	auto* p = range.data();
	for (uint32 i = 0; i < size; ++i) {
		p[i] = i;
	}
	return range;
}






CUDA_CALLABLE Inline vec3 vec_平方(const vec3& a) {
	return a * a;
};

CUDA_CALLABLE Inline float32 vec_和(const vec3& a) {
	return a.x + a.y + a.z;
};

CUDA_CALLABLE Inline float32 vec_平方范数(vec3 v) {
	v = vec_平方(v);
	return vec_和(v);
}





// 光线与圆球求交
CUDA_CALLABLE Inline float32 f_射线与点相交(const vec3& rayO, const vec3& rayD, const S_VecSphere& sphere) {
	vec3 oc = rayO - sphere.pos; 

	float32 a = vec_dot(oc, oc);
	float32 b = 2.0 * vec_dot(oc, rayD);
	
	float32 c = a - sphere.radius * sphere.radius;
	float32 h = b * b - 4.0 * c;
	if (h < 0.0) {
		return -1.0;
	}
	float32 t = (-b - sqrt(h)) / (2.0);

	return t;
}
CUDA_CALLABLE Inline float32 f_双向射线与点相交(const vec3& rayO, const vec3& rayD, const S_VecSphere& sphere) {
	vec3 oc = rayO - sphere.pos;

	float32 a = vec_dot(oc, oc);
	float32 c = a - sphere.radius * sphere.radius;

	float32 b = 2.0 * vec_dot(oc, rayD);

	float32 h = b * b - 4.0 * c;
	return (-b - sqrt(fabs(h))) / (2.0);

	if (h < 0.0) {

		b = 2.0 * vec_dot(oc, vec_反向(rayD));
		h = b * b - 4.0 * c;

		if (h < 0.0) return -1.0;
	}
	
	a = (-b - sqrt(h)) / (2.0);
	return a;
}


CUDA_CALLABLE Inline float32 f_直线到点的距离(const vec3& rayO, const vec3& 单位向量, const S_VecSphere& sphere) {
	vec3 oc = sphere.pos - rayO;
	float32 斜边长度 = vec_len(oc);
	
	float64 直角边长度 = 斜边长度 * 斜边长度 - sphere.radius * sphere.radius;
	直角边长度 = sqrt(abs(直角边长度));

	vec3 射线点 = 单位向量 * float32(直角边长度) + rayO;
	float32 检测距离A = vec_len(射线点 - sphere.pos);

	float32 距离A = sphere.radius - 检测距离A;
	if (距离A < 0) {
		射线点 = rayO - 单位向量 * 直角边长度;
		检测距离A = vec_len(射线点 - sphere.pos);

		return sphere.radius - 检测距离A;
	}
	return 距离A;




	float32 a = vec_dot(oc, oc);
	float32 b = 2.0 * vec_dot(oc, 单位向量);

	float32 c = a - sphere.radius * sphere.radius;
	float32 h = b * b - 4.0 * c;
	
	float32 t = (-b - sqrt(fabs(h))) / (2.0);
	return t;
}



//vec3 f_graph_距离最近坐标(vec3 原点, vec3* 坐标集, uint32 num);
Inline vec3 f_graph_距离最近坐标(vec3 原点, vec3* 坐标集, uint32 num) {
	float32 dis = vec_len(原点 - 坐标集[0]);
	vec3 最近点 = 坐标集[0];
	for (uint32 i = 1; i < num; ++i) {
		float32 dis_temp = vec_len(原点 - 坐标集[i]);
		if (dis_temp < dis) {
			dis = dis_temp;
			最近点 = 坐标集[i];
		}
	}
	return 最近点;
}




//取三角形AB边
Inline vec3 geometry_getTriangleEdgeAB(const S_三角形& triangle) {
	vec3 vecAB;
	vecAB.x = triangle.pointB.x - triangle.pointA.x;
	vecAB.y = triangle.pointB.y - triangle.pointA.y;
	vecAB.z = triangle.pointB.z - triangle.pointA.z;
	return vecAB;
}

//取三角形BC边
Inline vec3 geometry_getTriangleEdgeBC(const S_三角形& triangle) {
	vec3 vecBC;
	vecBC.x = triangle.pointC.x - triangle.pointB.x;
	vecBC.y = triangle.pointC.y - triangle.pointB.y;
	vecBC.z = triangle.pointC.z - triangle.pointB.z;
	return vecBC;
}

//取三角形AC边
Inline vec3 geometry_getTriangleEdgeAC(const S_三角形& triangle) {
	vec3 vecAC;
	vecAC.x = triangle.pointC.x - triangle.pointA.x;
	vecAC.y = triangle.pointC.y - triangle.pointA.y;
	vecAC.z = triangle.pointC.z - triangle.pointA.z;
	return vecAC;
}


Inline float f_graph_点到三角形平面距离(const S_三角形& triangle, const vec3& OutSidePoint) {
	vec3 vecAB = geometry_getTriangleEdgeAB(triangle);
	vec3 vecAC = geometry_getTriangleEdgeAC(triangle);

	vec3 vecPointToPlane =
	{ OutSidePoint.x - triangle.pointA.x, OutSidePoint.y - triangle.pointA.y, OutSidePoint.z - triangle.pointA.z };

	vecAB = vec_normalize2(vecAB);
	vecAC = vec_normalize2(vecAC);
	//计算平面法线
	vec3 PNormal = vec_cross(vecAB, vecAC);
	float dis = vec_normalize(&PNormal) + DEF_浮点精度补充;
	
	//得到点和平面的距离
	return vec_dot(vecPointToPlane, PNormal) / dis;
}

Inline float32 f_graph_点到平面距离(const vec3& point, const vec4& plane) {
	return plane.x * point.x + plane.y * point.y + plane.z * point.z + plane.w;
}

Inline bool f_graph_射线相交平面(const vec3& 起点, const vec3& 方向, const vec4& plane, float32& 相交距离) {
	float32 dist = vec_dot(plane, 方向);

	if (dist > 0.0f && dist < 1.0f) {
		return false;
	}
	else {
		相交距离 = -vec_dot(plane, 起点) / dist;
	}
	return true;
}

Inline bool f_graph_射线相交平面(const vec3& 起点, const vec3& 终点, const vec3& plane, vec3& out) {
	vec3 u = (终点 - 起点);

	float32 dist = -vec_dot(plane, 起点) / vec_dot(plane, u);

	if (dist > 0.0f && dist < 1.0f) {
		out = (起点 + u * dist);
		return true;
	}
	else return false;
}

Inline bool f_graph_射线相交平面2(const vec3& 起点, const vec3& 终点, const vec4& plane, vec3& out) {
	vec3 u = (终点 - 起点);

	float32 dist = -vec_dot(plane, 起点 + vec3{ plane.x, plane.y, plane.z } * plane.w) / vec_dot(plane, u);

	if (dist > 0.0f && dist < 1.0f) {
		out = (起点 + u * dist);
		return true;
	}
	else return false;
}

Inline float64 f_graph_两向量夹角弧度(const vec3& vec0, const vec3& vec1) {
	float64 角A = abs(vec_dot(vec0, vec1));
	float64 角B = vec_len(vec0) * vec_len(vec1);
	if(角B == 0) 角B = DEF_浮点精度补充;
	//return acos(角A / 角B);
	return acos( DEF_Min(DEF_Max(角A / 角B, 0), 1) );
}

float32 f_graph_两向量夹角弧度_export(const vec3& vec0, const vec3& vec1);

Inline float64 f_graph_两向量夹角弧度180(const vec3& vec0, const vec3& vec1) {
	float32 角A = vec_dot(vec0, vec1);
	float32 角B = vec_len(vec0) * vec_len(vec1);

	角B = acos(fabs(角A) / (角B + DEF_浮点精度补充));
	if (角A < 0) {
		角B = 3.14159265 - 角B;
	}
	return 角B;
}

Inline float64 f_graph_两向量夹角弧度(const vec2& vec0, const vec2& vec1) {
	float32 角A = (vec_dot(vec0, vec1));
	float32 角B = vec_len(vec0) * vec_len(vec1);
	return acos( DEF_Min(DEF_Max(角A / (角B + DEF_浮点精度补充), -1), 1) );
}
Inline float64 f_graph_两向量夹角弧度180(const vec2& vec0, const vec2& vec1) {
	float32 角A = (vec_dot(vec0, vec1));
	float32 角B = vec_len(vec0) * vec_len(vec1);
	角B = acos( DEF_Min(fabs(角A) / (角B + DEF_浮点精度补充), 1) );
	if (角A < 0) {
		角B = 3.14159265 - 角B;
	}
	return 角B;
}



Inline float64 f_graph_三角函数_临边计算斜边长度(const float64 临边B长度, const float64 角弧度) {
	return 临边B长度 / (cos(角弧度) + DEF_浮点精度补充);
}
Inline float64 f_graph_三角函数_斜边计算对边长度(const float64 斜边C长度, const float64 角弧度) {
	return sin(角弧度) * 斜边C长度;
}
Inline float64 f_graph_三角函数_斜边计算临边长度(const float64 斜边C长度, const float64 角弧度) {
	return cos(角弧度) * 斜边C长度;
}
Inline float64 f_graph_三角函数_对边计算斜边长度(const float64 对边A长度, const float64 角弧度) {
	return 对边A长度 / (sin(角弧度) + DEF_浮点精度补充);
}
Inline float64 f_graph_三角函数_对边计算临边长度(const float64 对边A长度, const float64 角弧度) {
	return 对边A长度 / (tan(角弧度) + DEF_浮点精度补充);
}
Inline float64 f_graph_三角函数_临边计算对边长度(const float64 对边A长度, const float64 角弧度) {
	return 对边A长度 * tan(角弧度);
}


float32 f_graph_两线段距离(const vec3& vec0, const vec3& vec1, const vec3& line起点, const vec3& line方向单位向量);
float32 f_graph_两线段距离(const vec3& vec0, const vec3& vec1, const vec3& line起点, const vec3& line方向单位向量, vec3& 相交点, float32 线段补偿距离 = 0);

Inline float32 f_graph_点到直线距离(const vec3& s, const vec3& dir, const vec3& point) {
	vec3 x = point - s;
	//float32 b = 2.0 * dot(x, dir);
	//float32 c = dot(oc, oc) - sphere.radius * sphere.radius;
	float32 l = vec_normalize(&x);
	//return f_graph_三角函数_斜边计算对边长度(l, vec_dot(x, dir));
	return f_graph_三角函数_斜边计算对边长度(l, f_graph_两向量夹角弧度(x, dir));
}


Inline vec3 f_graph_计算点到直线垂足坐标(const vec3& s, const vec3& dir, const vec3& point) {
	vec3 x = point - s;
	float32 l = vec_normalize(&x);
	float32 临边长度 = f_graph_三角函数_斜边计算临边长度(l, f_graph_两向量夹角弧度(x, dir));
	return s + dir* 临边长度;
}



vec3 f_graph_从向量计算旋转角(const vec3& 向量1, const vec3& 向量2);






// 反射光线
/*void reflectRay(inout vec3 rayD, in vec3 mormal)
{
	rayD = rayD + 2.0 * -dot(mormal, rayD) * mormal;
}

// 光线与圆球求交
float sphereIntersect(in vec3 rayO, in vec3 rayD, in Sphere sphere)
{
	vec3 oc = rayO - sphere.pos;
	float b = 2.0 * dot(oc, rayD);
	float c = dot(oc, oc) - sphere.radius * sphere.radius;
	float h = b * b - 4.0 * c;
	if (h < 0.0)
	{
		return -1.0;
	}
	float t = (-b - sqrt(h)) / 2.0;

	return t;
}

vec3 sphereNormal(in vec3 pos, in Sphere sphere)
{
	return (pos - sphere.pos) / sphere.radius;
}

// 光线与平面求交
float planeIntersect(vec3 rayO, vec3 rayD, Plane plane)
{
	float d = dot(rayD, plane.normal);

	if (d == 0.0)
		return 0.0;

	float t = -(plane.distance + dot(rayO, plane.normal)) / d;

	if (t < 0.0)
		return 0.0;

	return t;
}*/


/*Inline float f_点到平面距离(const S_三角形& triangle, const float OutSidePoint[3], float PNormal[3]) {
	float vecAB[3] = { 0.0,0.0,0.0 };
	float vecAC[3] = { 0.0,0.0,0.0 };
	geometry_getTriangleEdgeAB(triangle, vecAB);
	geometry_getTriangleEdgeAC(triangle, vecAC);
	float vecPointToPlane[3] =
	{ OutSidePoint[0] - triangle.pointA.x, OutSidePoint[1] - triangle.pointA.y, OutSidePoint[2] - triangle.pointA.z };

	vec_归一化(vecAB);
	vec_归一化(vecAC);
	//计算平面法线
	vec_叉积(vecAB, vecAC, PNormal);
	float dis = vec_归一化(PNormal);

	//得到点和平面的距离
	return vec_dot(vecPointToPlane, PNormal) / dis;
}*/

//Inline vec3 f_graph_从向量计算欧拉角(const vec3& dir) {
//	vec2 zdir = { dir.x, dir.y };
//	float32 z = f_graph_两向量夹角弧度180(zdir, { 1,0 });
//}





